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Phase diagram of supercooled water confined to hydrophilic nanopores 
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We present a phase diagram for water confined to cylindrical silica nanopores in terms of pressure, temperature 
and pore radius. The confining cylindrical wall is hydrophilic and disordered, which has a destabilizing effect 
on ordered water structure. The phase diagram for this class of systems is derived from general arguments, 
with parameters taken from experimental observations and computer simulations and with assumptions tested 
by computer simulation. Phase space divides into three regions: a single liquid, a crystal-like solid, and glass. 
For large pores, radii exceeding 1 nm, water exhibits liquid and crystal-like behaviors, with abrupt crossovers 
between these regimes. For small pore radii, crystal-like behavior is unstable and water remains amorphous 
for all non-zero temperatures. At low enough temperatures, these states are glasses. Several experimental 
results for supercooled water can be understood in terms of the phase diagram we present. 



Instances of water confined to nanoscopic dimensions 
are ubiquitous in nature and technology. For example, 
water confined to silica nanopores coated with catalyst 
is an efficient system for evolving oxygen, a first step to- 
wards artificial photosynthesis. 1 Another example, aqua- 
porin pores in biological membranes confine water to 
channels in such a way to control water content in cells P 
Further, silica nanopores are also used to inhibit freezing 
of water to enable exploration 3 5 of water's behavior at 
conditions where the bulk material would spontaneously 
crystalize. The extent to which behaviors of water con- 
fined in this way - to a long hydrophilic nanopore - re- 
flects behaviors of bulk water has been unknown. Here, 
we address this deficiency by using general theoretical ar- 
guments coupled with molecular simulation to construct 
the phase diagram for water at standard and supercooled 




conditions as a function of temperature, T, pressure, p, 
and pore radius. 

The class of systems we consider is illustrated in Fig. 
[T] which shows snapshots taken from a molecular simula- 
tion of water confined to a silica nanopore, details about 
which are given later. This confining pore is long and nar- 
row. Its walls are hydrophilic, but with atoms that are in 
a disordered arrangement, much like a typical arrange- 
ment of oxygen atoms in liquid water, except the atoms 
making up the pore are frozen in place and have a slightly 
larger space filling size than that of water. This static 
surface disorder inhibits crystal-like structures so that 
water adjacent to the surface is typically disordered too. 
We will show that the thickness of that adjacent layer is 
A ~ 2.5 A, about the diameter of one water molecule. 

By considering pore radii R p = R + A that are twice 
2.5 A or larger, there can remain a significant amount of 
confined water that is not part of the adjacent mono- 
layer. This interior water takes on ordered or disor- 
dered arrangements, depending upon temperature, pres- 
sure and pore radus. For small enough radii, the desta- 
bilizing influence of the amorphous layer causes interior 
ordered water to be unstable, even at temperatures far 
below standard freezing temperatures. For larger radii, 
where ordered states are thermodynamically stable, time 
scales over which an ordered structure might emerge can 
be very long, so long that the interior water may be- 
come glass. These are the features considered in this 
paper: bulk thermodynamic stability, competing inter- 
facial energetics and time scales to reorganize molecular 
structures. 



FIG. 1. Oxygen atoms of water (the red particles) confined 
in a molecular simulation to a nanopore of length L with a 
disordered hydrophilic surface chosen to mimic silica in the 
material named MCM-41.^ The radius of the pore is R + A, 
where A is the thickness of an amorphous water mono- layer 
adjacent to the pore surface. See text. 
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PHASE DIAGRAM 

Melting in a bulk macroscopic system coincides with 
a singularity in a free energy function. In a bounded 
system, like those we consider here, the transition is 
smoothed or altogether removed. Two relevant length 
scales associated with this behavior emerge from the mi- 
croscopic theory presented in the next section. The first, 
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FIG. 2. Phase diagrams for supercooled confined water, 
a) Melting and glass transition temperatures, T m (p,R) and 
T g (p, R) , respectively, as a functions of pressure p and cylin- 
der radius, R. b) Phase diagram at the constant pressure 
p — 1 bar. Triangular markers indicate melting temperatures 
measured experimentally. 7 Circles indicate melting temper- 
atures determined through our computer simulation, where 
errorbars indicate our uncertainty in T m (p, R) . Squares indi- 
cate glass transition temperatures measured experimentally 
with error estimates for R (not shown in figure) of about zb 4 
AP c) Phase dia gram for two different fixed radii R = 5.0 A 
(solid line) and R — 9.5 A (dashed line). Circles indicate an 
onset of thermal hysteresis in experimental density measure- 
ments with R « 5.0 A. 4 Error bars indicate our measure of 
uncertainty of where hysteresis begins. The square indicates 
an estimate of the calorimetric glass transition for a pore of 
approximately the same diameter.^ The error estimate stated 
in Ref. [8] is smaller than the size of the symbol. 



reflects competition between bulk energetics favoring or- 
der and interfacial energetics opposing order. Here, 7 
is the surface tension between the ordered crystal and 
the disordered liquid, and Ah is the heat of fusion per 
unit volume. The value, £ m « 0.21 nm for water, follows 
from the values of su rface t ension and heat of fusion for 
water-ice coexistence! 9 * 1 ^ 1 ^ Both 7 and Ah are pressure 
dependent, but the ratio Ah is pressure-independent 
to a good approximation.^ More is said about this fact 
later. 

When the radius R is significantly larger than £ m: a 
melting temperature remains finite. This temperature, 
T m (p,R), is defined as that where the free energy of an 
ordered structure equals the free energy of a liquid. Ac- 
cording to macroscopic thermodynamics, T m (p, R) fol- 
lows a Gibbs-Thompson equation (like the Kelvin equa- 
tion in the context of capillary condensation)!^ Specif- 
ically, T m (p,R) w T m (p)(l - i m /R% where T m (p) is the 
bulk melting temperature. This approximation describ- 
ing the reduction in melting temperature with increasing 
1/R is correct to the extent that £ m /R <C 1. The melting 
curve for small 1/R shown in Fig. [2] follows this equation. 

The second relevant length emerging from the theory 
manifests fluctuations that destabilize order. Specifically, 
fluctuations renormalize the first length to yield 



4= W(l-T s /T m ) ^0.91nm, 



(2) 



where T s stands for the temperature below which a bulk 
amorphous phase of water is unstable. It is gener- 
ally pressure dependent, but according to our simulation 
studies of one water mo del the ratio T s (p)/T m (p) is 
pressure independent. We therefore omit explicit refer- 
ence to its pressure dependence in Eq. [2] Experimen- 
tally, it is difficult to measure T s , so in order to es- 
timate it for water we rewrite the ratio as T s /T m = 
(T s /T pmax )(T pmax /T m ) where T pmax is the temperature 
of maximum density at low pressure (Ref. [14| uses the 
symbol T Q for that temperature). We write these ra- 
tios because we have found previously that T pmax rep- 
resents the relevant energy scale for supercooled water 
thermodynamics. Therefore, we expect that for any rea- 
sonable model of water T s /T pmax will be independent of 
the specific choice of model. As such, it can be extracted 
from simulation, with which we find it to be T s /T pmax = 
0.76.^ The second term, T pmax /T m , is a model depen- 
dent constant, often close to unity and its value for water 
is known experimentally to be l.OlP We use that value. 
Therefore, we predict that T s (latm) = 210 K for water. 
This prediction of a lower temperature limit to liquid 
stability is consistent with experimental observations of 
rapid spontaneous crystallization of water at 220KP^In 
addition, it yields the value 4 = 0.91 nm cited above. 

The Gibbs-Thompson correction to the bulk melting 
line is accurate only when order parameter fluctuations 
can be neglected. These fluctuations become dominant 
as R approaches 4- Specifically, in the next section we 
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derive 



T m (p,R)^T m (p) 1 
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for R > R C1 where R c is the positive root of the right 
hand side of Eq. [5] and is approximately equal to 
For R < R C1 T m (p, R) = 0. This expression is graphed in 
Fig. [2] The fluctuation contribution produces the precip- 
itous end to the melting line near R « 1 nm. The com- 
parison of data points and lines in Fig. [2] shows that our 
predicted behavior of T m (p, R) agrees well with observed 
calorimetry results for an order-disorder transformation 
of water in silica pores. 7 Equation [3] also agrees well with 
our molecular simulation results discussed later in this 
paper. 

The second surface shown in Fig. 2j T g (p,R), is de- 
fined to be the temperature below which the structural 
relaxation time, r, of supercooled liquid water would be 
larger than 100 s. Nonequilibrium perturbations taking 
place on shorter time scales, such as cooling rates in the 
range of O.lK/min to lK/min or faster, would take the 
liquid out of equilibrium. It is in that sense that T g (p, R) 
is the glass transition temperature. To estimate its be- 
havior, we note that r(T) generally follows the parabolic 
form below an onset temperature* * 18 J i.e., 

log (t/t ) = J 2 (1/T - 1/T ) 2 , T < T . (4) 

We adopt this expression together with 100 s = 
r(T g ,p, R). The reference time, r Q , the onset tempera- 
ture, T Q , and the energy scale, J, are generally functions 
of p and R. These functions can be determined from 
simulation and experiment. 

One such determination is that r Q for water is close to 
1 ps throughout the range of p and R we find relevant. 
Accordingly, Fig. [2] graphs 



T g (p,R)nT (p,R)/ [l + VUT Q (p,R)/J(p,R)\ . (5) 

From experiment and simulation of water, we can de- 
termine functional forms for J(p,R) and T (p,R). See 
Methods section. While each depends monotonically on p 
and i?, their systematic trends lead to the non-monotonic 
behavior T g (p,R) illustrated in Fig. [2] Most notable is 
how the slope of T g (p, R) with respect to p changes from 
small and negative when R is large to relatively large and 
positive when R is small. This variation in slope explains 
a few critical observations. 

In particular, recent calorimetry experiments prob- 
ing glassy relaxation in confined systems have estimated 
T g (p,R) for several pore sizepl at low pressures. The 
results of those observations (the squares in Fig. [5| coin- 
cide closely with our predictions for this glass transition 
temperature. Furthermore, Zhang et alW have, in ef- 
fect, located the glass transition temperature at higher 
pressures through their observation of hysteretic behav- 
ior for density in nano-pores upon cooling at a rate of 0.2 
K/min. Hysteresis occurs only because the system falls 



out of equilibrium. Data points from Ref. 4 (the circles 
in Fig. |2| fall close to our predicted glass transition tem- 
perature line. That reference attributes the hysteresis 
to something other than a glass transition, namely a hy- 
pothesized liquid-liquid transition.^ Previous work by us 
casts doubt on that possibility,^ leaving the glass transi- 
tion as a plausible explanation for pore sizes as small as 
those reported in Ref. |4j If the pore sizes were a fac- 
tor of 2 larger than estimated by those authors, our phase 
diagram indicates that hysteresis could also reflect time 
scales for nucleating an ordered crystal-like material. 

At temperatures below T g (p, i?), water may exist in 
more than one amorphous solid state. Preparations 
and transitions between these amorphous states are ir- 
reversible and therefore beyond the scope of this paper. 



DERIVATION OF PHASE DIAGRAM 

Our approach for analyzing remnants of first-order 
phase transitions in bounded systems begins by choos- 
ing a general phenomenological hamiltonian for an order- 
parameter field parameterized with experimental data. 
We then perform statistical mechanical calculations for 
the bounded systems based upon that hamiltonian, and 
we test assumptions in our analysis with atomistic simu- 
lations. Our strategy for examining the possibility of out- 
of-equilib rium transitions to glassy sates is to use scaling 
principle d 17 1 22 1 to bootstrap to the glass transition from 
knowledge of structural relaxation times at moderately 
supercooled conditions, and to use molecular simulation 
to test assumptions underlying that approach. 



Equilibrium 

We consider an energy functional or hamiltonian of an 
order parameter distinguishing a liquid-like state from a 
crystal-like state. For bulk water, the two states can be 
distinguished with a global order parameter like Stein- 
hardt, Nelson and Ronchetti's Qq variable. 23 Complex 
fields for local order parameters could be used too. Bro- 
ken symmetry for either Qq or a phase of a complex field 
does not occur for confined systems like those we con- 
sider here. Therefore, we choose to distinguish liquid-like 
states from more ordered crystal-like states in terms of a 
local order field that is real. There are many such mea- 
sures suitable for our purpose. As a specific example, our 
choice of order parameter could be 
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Here, the sum over j G nn(i) includes only the 4 near- 
est neighbor oxygens of the ith oxygen, and l6m(0ij? Oij) 
is the i = 6,m spherical harmonic function associated 
with the angular coordinates of the vector — rj joining 
molecules i and j measured with respect to an arbitrary 
external frame. This particular order parameter is large 
in proportion to the concentration of water molecules 
with neighbors having the same orientations of neighbor- 
ing bonds as does the molecule itself. The quantity qu q is 
its non-zero value for the bulk liquid. Past experience has 
shown that using the £ = 6 spherical harmonics with 4 
nearest neighbors is particularly useful for detailing local 
structure in water ™ 

With this or some similar order-parameter field, we 
choose the energy to have the following form 
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wq + uq 
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where &b is Boltzmann's constant, q(r) is the deviation 
of the order parameter field from its uniform value for 
bulk liquid water, and a = ao(T — T s ), with T s being the 
temperature below which the bulk liquid is unstable. The 
parameters ao, w, u, m and T s are positive constants that 
depend upon pressure but are independent of tempera- 
ture. They can be determined in terms of the measured 
properties of bulk water, all as specified below. The zero 
of energy is that of the disordered amorphous material. 

The most significant feature of this phenomenological 
energy functional is the presence of only one field, specifi- 
cally one that refers to local order, and that local molecu- 
lar density is explicitly absent. We adopt this feature for 
two reasons. First, our prior simulation work indicates 
that the reversible free energy function for condensed wa- 
ter has no more than one amorphous-state basin and no 
more than one crystal-state basin at the conditions we 
consider.^ Second, we show below that we do not need 
to invoke the possibility of two liquid states to explain 
the experimental data we set out to interpret. In our 
picture, therefore, molecular density introduces no addi- 
tional phase-transition-like behavior. To the extent our 
picture is accurate, integrating out density therefore af- 
fects only the values of the parameters already included. 

Two other significant features of the energy functional 
is the truncation at fourth order in the order parame- 
ter and the neglect of inhomogeneity beyond the square 
gradient term. The first feature limits our treatment to 
no more than two distinct reversible phases. As noted 
above, we believe this limitation is acceptable for the 
conditions we consider. The second feature limits our 
treatment to a long wavelength description of interfaces. 



It is the simplest description for estimating the role of 
surface energetics. 24 In our simulations for the crystal- 
like phases, we find instantaneous domains with faceted 
surfaces. These features are smoothed by averaging over 
the disorder imposed by the pore wall. The energetics 
we aim to describe with the square-gradient term is for 
interfaces averaged over the full length of the pore. The 
Methods section presents tests of this idea. 

Given the energy functional, the free energy, F(T,_p), 
is determined by the partition function, 



F(T,p) 
k B T 



= - In j Vq(r) exp {-H[q(r)]/k B T} , (11) 



where the dependence upon pressure, p, enters through 
the parameters in the model, as discussed above and de- 
tailed further below. 



Mean field treatment 

Mean field theory identifies the mean value (q(v)) as 
the function q(r) that minimizes H[q(r)] subject to the 
boundary conditions imposed by the cylindrical pore. 
Specifically, the free energy in this approximation is 



Fmf = U[(q(r))] , 



where 



6H 
,%(r)> 



0. 



(12) 



(13) 



For boundary conditions applied to solving Eq. [13j we 
assume the effects of the pore are two-fold. First, we 
assume the pore confines water to a cylinder of radius 
R p . Second, we assume that disorder of the pore's con- 
fining hydrophilic surface induces liquid-like behavior in 
adjacent water, making 



(q(r)) = 0, for |r| > R p - A = R, 



(14) 



where A is the thickness of the amorphous boundary 
layer. Simulation results for various temperatures and 
pore radii (see Methods section) show that crystal-like 
domains are limited to an inner cylinder of radius R = 
R p — 2. 5 A where R p is the mean distance from the center 
of the pore to the silica wall. Therefore, we take A ~ 2. 5 A 
for all temperatures and pore radii. 

In this picture, surface energetics controlling the non- 
trivial behavior of (q(r)) is determined by the interface 
between liquid water and ice. The silica- water interac- 
tions are irrelevant except for producing a layer of dis- 
ordered water and therefore imposing the boundary con- 
dition Eq. 14 One important consequence is that the 



parameter m is determined by the interfacial energy of 
ice in contact with liquid water. 

For an analytical solution to Eq. [l3j we consider R to 
be very large. To leading order in 1/i?, the solution to 
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Eq. 13 is that of a one dimensional interface that yields 
a mean field free energy per unit volume 



with 
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where 



dF M F/dq = 0, 



(15) 



(16) 



with q being the mean-field estimate of (g(r)) for r at 
the center of the cylinder. With numerical solutions, we 
have checked that terms beyond linear in 1/R would be 
significant contributors to Fmf(q)/V only for R so small 
that fluctuation corrections to the mean field approxima- 
tion are dominant (see below). At that stage, corrections 
to Eq. 15 due to growing curvature are irrelevant. 

At a point where a disordered state, q = 0, coexists 
with an ordered state, q = q xt \ > 0, the mean field con- 
dition for coexistence is 



^Mf(O) = i^MF^xtl) : 



(17) 



which has a solution at a temperature Tmf(p, R), and 
entropy difference per unit volume given by 
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d^MF(gxtl) 

dT 



(18) 



T=T MF (p,R) 



In the limit R — >• oo, this coexistence should coincide 
with bulk freezing transition because it is a first-order 
phase transition where fluctuation effects are not impor- 
tant. Accordingly, we associate Tmf(p,R oo) with 
the bulk freezing temperature, T m (p), and As(p,R — >• 
oo) = As(p) with the entropy change between water 
and ice, i.e., As(p) = —Ah/T m (p). These connections 
to the bulk melting transition together with the corre- 
sponding mean field approximation for surface tensionf" 
7 = J^ 1 [2mf(q)]~ 1/2 dq, allow us to identify all rele- 
vant combinations of parameters in the model in terms 
of experimentally observed properties of bulk water!^ 
Specifically, after some algebra Eqs. [I] and 16 yield a 
mean-field expression for the melting surface, 



25 



T MF (p,R)=T m (p)(l-e m /R). 



(19) 



The mean-field approximation Tmf(p,R) ~ T m (p,R) is 
identical to the macroscopic Gibbs-Thompson estimate 
noted earlier. 



Role of fluctuations 

To estimate the effects of fluctuations, we evaluate 
AF(q) = F(q) — Fmf(q) in a Gaussian approximation. 
That is, 

AF(q) = -k B T In J Vq(r) exv{-AH[q(r)]/k B T} , 

(20) 



AH[q(r)] « ^ J^[5q(v)f^m\V5q(v)\ 2 } , (21) 

where n = a — 3wq + 6uq 2 . This approximation 
to AH[q(r)] comes from expanding T-L[q(r)] through 
quadratic order in Sq(r) = q(r) — q. 

The geometry of the system plays a role through the 



Laplacian in Eq. 21 For the cylinderical boundary con- 
ditions we consider, evaluation of the Gaussian integral 
prescribed by Eqs. 20 and 21 can be done by using zeroth 



order Bessel functions with the limits of integration re- 
stricted to allow fluctuations of wavelengths up to 2n/R. 

The resulting approximation to the free energy can be 
used to estimate the temperatures and pressures where 
the ordered and disordered materials have equal statisti- 
cal weight, i.e., where Eq. 17 is satisfied but with Fmf(q) 
replaced with the fluctuation corrected F(q). After some 
algebra, we find for R < R c 



T m (p,R) 
T m (p) 



1 



R 



(22) 



8tt(R - QR 



R 



where we have noted the order of neglected term. This 
term, 0{t 2 jR 2 ), refer to a curvature correction that 
would distinguish slab and cylinder geometries. For wa- 
ter at the conditions we consider, the dominant contribu- 
tion for small R is due to £ s / (R—£ s ) being large. As such, 
Fourier components rather than Bessel functions could 
have been used to diagonalize the determenent for the 
Gaussian integral, and equivalently, the partition func- 
tion we consider is dominated by its largest eigenvalue. 

The vanishing of crystal-like stability predicted in this 
way, where T m (p,R) — >> for R — >> R C1 is essentially 
a Ginzburg criterion.^ The length £ s is close to but 
necessarily smaller than this smallest radius, R C1 where 
crystal-like states can be stable. 



Nonequilibrium 

To evaluate the glass transition temperature from 
Eq. [4j we must determine r Q , T Q (p, R p ), and J(p,R p ). 
These parameters control very long-time relaxation, but 
they can be accessed through computation and exper- 
iment that measure relatively short time behavior.^ 
The Methods section describes our handling of exper- 
imental and simulation data to obtain these parame- 
ters. For bulk water, measured relaxation times yield 
To(latm) = T « 271 K, while the mW model used in 
our simulation yields T Q w 234 K; similarly, for bulk wa- 
ter J(l atm) = J ~ 7.5 T Q , while the mW model used in 
our simulation yields J « 6.3 T Q . 

In creating our phase diagram, we use the real-water 
values for these quantities. Nevertheless, the compari- 
son between these quantities for real water and for the 
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mW model give us confidence in using simulation to esti- 
mate quantities not available from experiment. In partic- 
ular, because liquid structure of mW water is virtually 
identical to that of real water, 32 we expect that rela- 
tive dependence upon R p can be accurately estimated 
with the simulation. The dependence we find in that 
way for R p > 5 A is T Q (R p ) w T G [1 + (6.0 A/i? p ) 2 ], and 
J(Rp) « J(l-4.4 k/Rp), where T Q (R P ) and J(# p ) stand 
for the low pressure values for T Q (p, i? p ) and J(p, i? p ), re- 
spectively. 

For the pressure dependence of these quantities, we 
rely on experimental measurements of relaxation times 
at R p ~ 7.5 A. 3 That data allows us to estimate first and 
second derivatives with respect to pressure, leading us to 
write 

J(p, Rp) « J(iip) + 490 (K/kbar 2 )^ 2 (23) 

and 

T (P,R P ) « r o (iip) - 26(K/kbar)p, (24) 

where J(R P ) and T Q (Rp) are given in the paragraph 
above. 



+ 4 + OOA##0 Ref.27 
Ref. 3 ♦ Ref. 28 ■ Ref. 29 mW model 
♦ AAA □ Ref. 30 ^ Ref. 31 I|IDAA 




J(p,R p )[l/T-l/T (p,R p )] 

FIG. 3. Collapse of confined liquid water relaxation times for 
different pore sizes and at different external pressures. Pri- 
mary graph is for T > T m (p,R). The data is from our sim- 
ulation results and experimental resultsP^^ Inset graph is 
for relaxation times of crystal-like state, i.e., T < T m (p,R). 
It includes our simulation results and remaining data taken 
from Refs[3l|27H3Ql 



We have checked that these algebraic forms accurately 
extrapolate from the low-pressure values of the mW 
model at finite R p to the high-pressure values for the 
mW model at 1/Rp -» 0. 

Using these forms for the transport parameters, and a 
value of t = 1 ps, we can collapse experimental data and 
our simulation results across pressures and pore sizes. 
Figure [3] shows this collapse where we have restricted 
the data to include only equilibrium liquid relaxation, 
i.e., T > Tm (v. R ) . Figure. [3] includes data from both 
experimentP^ZHHIl an d from our simulation study. While 
external pressures can be accurately controlled and re- 
ported, errors in pore sizes are large. 33 The Methods sec- 
tion discusses estimates of R p from exper imen tal data. 

Previous simulation^ and experimentP^H^ studies 
have indicated that confined liquid water undergoes an 
abrupt crossover in the temperature scaling of its relax- 
ation time. This crossover is a manifestation of a tran- 
sition between the liquid and crystal-like regimes, which 
we turn to now. 



RELAXATION OF CRYSTAL-LIKE STATES 

Molecular motion of crystal-like states in confinement 
takes place preferentially near the water pore interface, 
where the molecules are locally disordered. Like defect 
motion in a bulk crystal, though with a smaller barrier 
due to the presence of the interface, the temperature de- 
pendence of such motion is expected to be Arrhenious. 
The inset of Fig. [3] shows experimental and simulation 
data for the average time for a particle to displace one 
molecular diameter as a function of temperature at con- 
ditions where water in the core is ordered. Plotted is ex- 
perimental data for T < T m (R) and simulation results for 
R = 17.5 A. We find that this motion is activated with 
a barrier of approximately 20 kJ/mol, and an attempt 
frequency l/r m ~2 ps -1 . There is negligible dependence 
on radius of confinement within the range considered. 

By combining the information of the phase diagram 
with our understanding of the mobility in each state 
we can predict the observed equilibrium behavior of the 
relaxation time. We find that there are three differ- 
ent pore size regimes, each with a distinct tempera- 
ture dependence of r. These regimes are highlighted 
in Fig. [4] First, for larger pores, R > 2£ s , the on- 
set to glassy dynamics is close to T m (p,R), therefore an 
equilibrium measurement should show little temperature 
dependence for T > T m (p,R) and an Arrhenious tem- 
perature dependence for T < T m (p,R) reflecting the re- 
laxation behavior of the crystal-like states. For smaller 
pores close to but larger than 4, the onset temperature is 
greater than T m (p, i?), therefore an equilibrium measure- 
ment should show parabolic temperature dependence for 
T > T m (p,R), and a crossover to Arrhenious behavior 
for T < T m (p, R). For very small pores, R < £ s (but still 
larger than a molecular diameter), T m (p, R) = there- 
fore an equilibrium measurement should show parabolic 
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temperature dependence for T < T Q . Figure [4] shows that 
each of these regimes are observed both in simulation and 
in experiment. 

We are not the first to suggest that the 
abrupt crossover in relaxation might be linked to 
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FIG. 4. Transport behavior for different pore sizes as indi- 
cated. Data are from simulation results (black points) and 
experimental results (red points) . Lines are predictions based 
on our phase diagram and scaling relations. Vertical dashed 
line in the middle panel locates the boundary between liquid 
and crystal- like states, i.e., where T = T m (p, R). 



crystallization.^ Some may have disregarded this possi- 
bility due to the absence of a freezing peak in the heat 
capacity, measured by differential scanning calorimetry. 
Our analysis shows that the absence of this peak is due 
to the pore size being close to £ s . When R « £ s , the 
ordering transition is smeared due to large structural 
fluctuations. As a consequence, there will be no sudden 
heat release. This explanation is consistent with a recent 
differential scanning calorimetry study that observed 
only partial crystallization for a pore size R p = 10.5, 
with the accompanying heat capacity peak being of the 
order of the magnitude of the maximum liquid state 
heat capacity. 8 

Experimentally determined vibrational density of 
states for confined supercooled water differs significantly 
from that of bulk ice, even at points in the phase diagram 
where we predict the presence of crystal-like behavior. 
This difference in density of states is expected because 
the domain of crystal-like behavior in the confined sys- 
tem is surrounded by a pre-melting layer, which in turn is 
surrounded by a layer of complete disorder. These layers, 
discussed in the Methods section, encompass a significant 
fraction of the total system, a fraction that grows with 
decreasing pore size. Further, even away from the dis- 
ordered pore wall, crystal-like behavior in confinement 
exhibits a high concentration of stacking faults, 36 which 
will further modify the density of states. 



METHODS 

Molecular simulation model 

The molecular dynamics simulations used to test our 
theoretical approximations and estimate the magnitudes 
of some differential changes employ the mW model of 
water Recently proposed by Molinero and Moore, 
this model has proven to yie ld a good description of 
water in the liquid phase J^"^, it reproduces many 
structural transformations seen in experiment and in 
other mo dels of water (including freezing into an ice-like 
structure) l 14 l 4Q t EH and it exhibits the characteristic ther- 
modynamic and dynamic anomalies of water (i.e., density 
maximum, heat capacity increase, diffusion maximum, 
and so forth) PEH 

To model the hydrophilic pores of MCM-41-S we have 
followed a procedure similar to that found in Ref. [43] 
The pore configurations are obtained by quenching a 
high temperature liquid configuration of silica. To create 
the cylindrical geometry, we extract from the simulation 
box all atoms whose centers lie within a circle of radius 
= l(%i-%c) 2 + (yi-y c ) 2 } 1/2 where (x c ,y c ) is the cen- 
ter of the simulation box and (xi,yi) is the coordinate 
vector for particle i. The remaining atoms are tethered 
to their initial conditions by a spherically symmetric har- 
monic potential with a spring constant, 50 kcal/mol A 2 . 
This procedure yields a mean surface roughness for the 
pore walls in good agreement with that estimated from 
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FIG. 5. Average value of the global orientational or- 
der parameter, (Qe), as a function of temperature for dif- 
ferent pore sizes. In the main figure, (Qq) is rescaled 
and plotted as a function of T — T*(R), where T* is 
the temperature at which | d(Qo)/dT | is maximal. In- 
set shows the same, but not rescaled. Different mark- 
ers correspond to different pore size systems, with R p = 
2.5 A + R = 20.0 A (blue Q) , 17.5 A (green □) , 15.0 A(cyan A) , 
12.5 A (red V), 10.0 A (grey O) and 7.5 A (black o). 



Determination of T m (R) from simulation 

In order to determine the low pressure melting tem- 
perature in confinement we calculate the temperature de- 
pendence of the orientational order parameter Qq defined 
as 



Qe 



6 N 

E E 

m= — 6 i,j 



1/2 



J * 



(25) 



following Ref. |23j Rather than a local measure of or- 
der, we use this global measure of crystallinity to be sure 
we are distinguishing crystal from liquid. Figure [5] shows 
the mean value of orientational order parameter, (Qq), as 
a function of temperature for six pore radii, i? p , ranging 
between 20.0 A and 7.5 A. For the range of temperatures 
we consider, pores with R p < 12.5 A never show pseudo 
long range order. Pores with R p > 12.5 A do show 
pseudo long range order. The presence of the amorphous 
interface ensures that (Qq) converges relatively quickly 
in comparison to the bulk where large nucleation bar- 
riers would separate the ordered and disorder states at 
coexistence. 

Apart from shifting the coexistence temperature, the 
pore radius changes the maximum value Qq(T) obtains 
in the ordered state. This is due to the increased weight 



experiment on MCM-41-S materials.^ We have consid- 
ered pore sizes in the range R = 5.0 A- 17.5 A. All 
pores are length L = 220 A to approach a regime where 
L > R. For initializing the combined water, pore system, 
water molecules are arranged in a hexagonal ice lattice 
at a density of 0.98 g/cm and placed within the pore 
with the crystallographic c-axis parallel to the length of 
the pore. 

All of the molecular dynamics trajectories were propa- 
gated using the LAMMPS package^ and a Nose- Hoover 
thermostat with constant number of particles, N, vol- 
ume, V, and temperature T. We fill approximately 90% 
of the length of the pore. Water organizes spontaneously 
with an interface separating the remaining 10% empty 
pore from the condensed phase (either liquid or crystal- 
like). With this procedure, we simulate the condensed 
material at a low pressure (effectively p « 0) in coexis- 
tence with vapor. 

We have adopted an interaction potential for a single 
site model of silica that is the same form as the mW 
model, but we have rescaled the interaction strength and 
particle diameter. Compared to the parameters used 
in the mW model, these are e s iii ca /e m w = 1-15 and 
Csiiica/^mW = 1-05. The increased interaction strength 
ensures a hydrophilic surface and the increased particle 
size frustrates local favored structures. 
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FIG. 6. Mean number density, (p(r)), and orientational order 
density, (q(r)), for water confined to cylindrical pores, r is 
the radial position of the pore system. Different color solid 
lines are our mW model simulation results for different pore 
size systems, with R p = 20.0 (blue), 17.5 (cyan), 15.0 (black) 
and 12.5 (green) A. All are computed at T « T m (p,R) and 
p « 0. The red dashed line is the prediction from the mean 
field theory, Eq.13, using parameters for 7i[q(r)] found with 
the mW model. A = 20 — R p is used to shift r to facilitate 
comparison of different radii pores. Notice that the order 
remains absent for density within A = 2.5 A of the pore wall. 
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that the amorphous boundary layer has on the volume in- 
tegral for decreasing R and fixed thickness A. The curves 
calculated for different R can be collapsed by multiply- 
ing Qq by a scaling factor. This factor accounts for the 
total volume available within pore compared to that for 
crystal-like states, (R+2X) 2 /R 2 , and subtractes the value 
for a disordered system of N water molecules. 



Disorder width, curvature and pre-melting layer 



Figure 6 shows the profiles for mean density, (p(r) 



and mean order parameter, (q(r)), as defined in Eq. [6] 
These curves are obtained from the mW model by sim- 
ulation, and from our theory by solving Eq. [l3]with the 
parameters appropriate for the mW model. The thermo- 
dynamic conditions considered are where the crystal-like 
state is stable. Several pore diameters were studied, and 
the illustrated results are typical. In each case, the simu- 
lation yields a disordered layer of non-zero particle den- 
sity and of thickness A = 2.5 A adjacent to the pore wall. 
This thickness of the disorder layer is in good agreement 



with the value inferred from fitting Eq. 19 to experimen- 
tal dataP 

The simulations also show oscillations in both the den- 
sity and the order parameter. These oscillations reflect 
the size of the particles in the simulated model. By con- 
struction, the square-gradient theory does not contain 
these oscillations. Nevertheless, the rise in the mean or- 
der parameter from its disordered value at the wall to its 
crystal-like value in the center of the pore is consistent 
with those of the simulation when coarse grained over 
a particle diameter. The general agreement of the pro- 
files calculated with our molecular dynamics simulations 
with those calculated neglecting curvature corrections in- 
dicates those corrections are small. 

The amplitude of the oscillations in the mean order 
parameter obtained from the simulation results are rela- 
tively small, typically 10% of the mean, except at the very 
center of the pore where statistics is unreliable. Away 
from the center, the oscillations are especially small in 
comparison to those that would be found in an ordered 
crystal. The amplitudes are diminished from those of a 
crystal due to the average over disorder along the length 
of the tube. 

The width of the interface exhibits a slight tempera- 
ture dependence. For larger pores, however, the situa- 
tion changes. As the radius grows beyond the conditions 
treated here, the coexistence temperatures will tend to- 
wards the bulk melting temperature. A pre-melting layer 
between the disordered surface and the crystal will then 
become large and strongly sensitive to temperature!^ For 
macroscopic systems, this pre-melting width diverges as 
T approaches the melting temperature. With the equa- 
tions we use in our theory, this behavior is isomorphic to a 
liquid- vapor wetting transition. It is is a general behavior 
accompanying any first-order transition with appropriate 
boundary conditions.^ 



The theory we have employed to describe pre-melting 
in a finite system would not only seem easily generaliz- 
able to treating pre-melting profiles of bulk ice in contact 
with its vapor or liquid, 46 it would also seem applicable 
to stability and thermodynamics of nano-clusters of ice,E3 
and to nucleation of ice on atmospheric aerosols.^ It 
might also be generalizable to describe ordering of water 
in cold micro emulsions like those recently considered by 
Tanaka and co-workers Indications of order-disorder 
phenomena occurring in the finite water-rich domains of 
those systems have been interpreted in terms of a doubt- 
ful liquid-liquid transition in supercooled water. Based 
upon what we have derived in this paper, we believe a 
more natural explanation of Tanaka' s observations will 
be found in terms ice-water equilibrium and the effects 
of confinement on that phase equilibrium 



Turnbull relation 

To test the applicability of Turnbull's j/Ah « 
constant, we have calculated the surface tension and en- 
thalpy of fusion as functions of temperature.^ See Fig.[7[ 
Here, Ah(T m ) = 5.4 kJ/mol and 7 (T m ) = 35.3 mJ/m 2 . 

To make that figure, we have determined the enthalpy 
of fusion by calculating the average enthalpy density dif- 
ferences at coexistence, (h)u q — (ft)xti- Similarly, we have 
determined the surface tension by calculating the free en- 
ergy as a function of Qq, using the umbrella sampling pro- 
cedure described in Ref.[14|for N = 216 particles at a con- 
stant pressure, p = 1 atm. See Fig. [7[ The surface ten- 
sion is then obtained by taking the difference between the 
free energy at the top of the barrier and at its stable coex- 
isting basins. Specifically, 7 = AF(Qq)/ L 2 where AF is 
the interfacial free energy calculated by first preforming 
a Maxwell construction to place the system at coexis- 
tence at the different temperatures, and L = (N '/ p) 2 / 3 P^ 
This procedure is exact in the limit of N — » 00. We 
have checked that we closely approach the limiting value 
by studying several system sizes up to N = 1000 parti- 
cles. This surface tension is an effective surface tension 
obtained by integrating over all distinct crystallographic 
faces and agrees well with that obtained from a recent 
nucleation study.™ 



T m (p, R) for water and the mW model 

For our calculations here we assume that the bulk 
melting line can be accuratly approximated by T m (p) = 
T m [l-pC + C(p 2 )}. The coefficient C is related to heat 
of fusion and the change in volume between water and ice 
determined at ambient pressure, as derived through the 
Clapeyron equation. For water C = 0.026 kbar -1 (Ref. 
|9j) and for the mW model C = 0.01 kbar" 1 P This dif- 
ference in slope between the mW model and real water 
is due to the mW model over estimating the density of 
Ice Ih. 32 However, by defining a pressure scale in units 
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FIG. 7. Validation of the Turnbull relation, j/Ah ~ const for the mW model, (a) Temperature dependence of the mean 
enthalpy in the bulk liquid and crystal states, (b) Liquid-crystal surface tension, 7, and the enthalpy of fusion, Ah, divided 
by their values at coexistence, T m = 274 K. (c) Free energy computed for a N = 216, p = 1 atm bulk system illustrating data 
with which coexistence and surface tension is determined. 



of C, the equation of state of water and the mW model 
can be related. 

The mW model we calculate £ m and £ s to be 2.40 A 
and 8.16 A respectively. These are slightly different than 
what we find for real water, and the differences account 
for the differences between our predicted melting line for 
real water and the calculated melting temperature of the 
mW model for 1/R = 0.1 A . See Fig. [2] The melting 
lengths for both the mW model and experiment agree 
with previo usly rep orted values based on fitting melting 
data to Eq.[l9f 
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Transport regression analysis 



We preformed a regression analysis on the algebraic 
forms used for J(p,R p ) and T Q (p, i? p ), Eqs. 



23 



and 



24 



See Fig. [8] The top two panels concern the dependence 
upon i? p , and the data for relaxation times is obtained 
from our molecular dynamics simulations of the confined 
mW model. The bottom two panels concern the pressure 
dependence, and the data for relaxation times is taken 
from experiments on confined water with R p = 7.5 AP 
The dashed lines are our algebraic fits, where the cor- 
relation coefficients indicate a certainty of 1% or better. 
Relaxation times for the molecular dynamics simulations 
were determined by calculating the mean time for a par- 
ticle to displace one diameter, 



FIG. 8. Onset temperature, T (p, R) and energy scale, 
J(p,R). The top two graphs are low pressure data found 
from simulations of the mW model. The bottom two graphs 
are high pressure data found from fitting experimental results 
of Ref. 3 In that case, the experiments report R p — 7.5A.The 
dashed lines are the curves obtained with Eqs. |23| and |24| 



Determination of R p from experimental data 

While nominal nanopore radii are routinely reported in 
the literature, it is difficult to obtain a reliable estimate 
of R p for R p < 1 nm. Different techniques yield a range 
of different sizes Most commonly, pore sizes are in- 
ferred from a Barrett- Joyner-Halenda (BJH) analysis.^ 
This method amounts to measuring a nitrogen absorp- 
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tion isotherm, and is thus an indirect measure of size. 
Mancinelli et al have demonstrated that this method can 
yield significant errors, up to 200 %. For example, using 
refined neutron scattering data and mass balance calcu- 
lations, Mancinelli et al. have estimated a likely range of 
pore sizes for the system studied by Ref. [3] to be between 
7.5A - 12.6AP while the BJH method yields 7.5A±2A. 

Using the bounds provided by Ref. 33 as reliable es- 
timates of possible errors, we find that we can collapse 
the experimental transport data, but such a collapse can- 
not be obtained within the errors reported from the BJH 
method. The pore sizes inferred from this collapse in- 
dicates that the BJH method systematically underesti- 
mates pore sizes Previous studies claiming to study the 
same pore sizes have observed widely different behavior. 
For instance, Ref. 3 report a pore radii R p = 7.5A, and 
measure a relaxation time that is never larger than 10's of 
nanoseconds. Reference |8] in one experiment also report 
using a pore of radii R p = 7. 5 A and measure thermal 
signatures indicative of a glass transition implying relax- 
ation time on the order of seconds. In light of our results 
detailing the different transport regimes that can occur 
for slightly different pore sizes, the implications of the er- 
rors associated with the reported values of the pore size 
become significant. 
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